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Abstract 

The concept of the limiting step is extended to the asymptotology of multiscale reaction networks. Complete theory 
for linear networks with well separated reaction rate constants is developed. We present algorithms for explicit 
approximations of eigenvalues and eigenvectors of kinetic matrix. Accuracy of estimates is proven. Performance of 
the algorithms is demonstrated on simple examples. Application of algorithms to nonlinear systems is discussed. 
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1. Introduction 

Most of mathematical models that really work 
are simplifications of the basic theoretical models 
and use in the backgrounds an assumption that 
some terms are big, and some other terms are small 
enough to neglect or almost neglect them. The closer 
consideration shows that such a simple separation 
on "small" and "big" terms should be used with 
precautions, and special culture was developed. The 
name "asymptotology" for this direction of science 
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was proposed by Kruskal (1963), but fundamental 
research in this direction are much older, and many 
fundamental approaches were developed by I. New- 
ton (Newton polyhedron, and many other things). 

Following Kruskal (1963), asymptotology is "the 
art of describing the behavior of a specified solution 
(or family of solutions) of a system in a limiting case. 
... The art of asymptotology lies partly in choosing 
fruitful limiting cases to examine ... The scientific el- 
ement in asymptotology resides in the nonarbitrari- 
ness of the asymptotic behavior and of its descrip- 
tion, once the limiting case has been decided upon." 

Asymptotic behavior of rational functions of sev- 
eral positive variables ki > gives us a toy-example. 
Let 

R(k x , ...k n ) = P(ki, . . . k n )/Q(ki, ...k n ) 
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be such a function and P, Q be polynomials. To de- 
rive fruitful limiting cases we consider logarithmic 
straight lines In k{ = 0i£, and study asymptotical be- 
havior of R for £ — >• oo. In this asymptotics, for al- 
most every vector (Oi) (outside several hyperplanes) 
there exists such a dominant monomial R ao (k) = 
A Yl t k T that R = R oo + o(Roo). The function that 
associates a monomial with vector (Oi) is piecewise 
constant: it is constant inside some polyhedral cones. 

Implicit functions given by equations which de- 
pend on parameters provide plenty of more inter- 
esting examples, especially in the case when the 
implicit function theorem is not applicable. Some 
analytical examples are presented by Andrianov & 
Manevitch (2002) and White (2006). Introduction 
of algebraic backgrounds and special software is pro- 
vided by Greuel & Pfister (2002). 

For a difficult problem, analysis of eigenvalues 
and eigenvectors of non-symmetric matrices, Vishik 
& Ljusternik (1960) studied asymptotic behavior of 
spectra and spectral projectors along the logarith- 
mic straight lines in the space of matrices. This anal- 
ysis was continued by Lidskii (1965). 

We study networks of linear reactions. For a linear 
system with reaction rate constants fcj all the dy- 
namical information is contained in eigenvalues and 
eigenvectors of the kinetic matrix or, more precisely, 
in its transformation to the Jordan normal form. It is 
computationally expensive task to find this transfor- 
mation for a non-symmetric matrix which is usually 
stiff (Golub & Van Loan (1996)). Moreover, the an- 
swer could be very sensitive to the errors in constants 
fcj. Nevertheless, it appears that stiffness can help 
us to find a robust approximation, and in the limit 
when all constants are very different (well-separated 
constants) the asymptotical behavior of eigenvalues 
and eigenvectors follow simple explicit expressions. 
Analysis of this asymptotics is our main goal. 

In our approach, we study asymptotic behavior 
of eigenvalues and eigenvectors of kinetic matrices 
along logarithmic straight lines, In fcj = in the 
space of constants. We significantly use the graph 
representation of chemical reaction networks and 
demonstrate, that for almost every vector (Oi) there 
exists a simple reaction network which describes 
the dominant term of this asymptotic. Following 
the asymptotology terminology (White (2006)), we 
call this simple network the dominant system. For 
these dominant system there are explicit formulas 
for eigenvalues and eigenvectors. The topology of 
dominant systems is rather simple: they are acyclic 
networks without branching. This allows us to con- 



struct the explicit asymptotics of eigenvectors and 
eigenvalues. All algorithms are represented topolog- 
ically by transformation of the graph of reaction (la- 
beled by reaction rate constants) . The reaction rate 
constants for dominant systems may not coincide 
with constant of original network. In general, they 
are monomials of the original constants. 

This result fully supports the observation by 
Kruskal (1963): "And the answer quite generally 
has the form of a new system (well posed prob- 
lem) for the solution to satisfy, although this is 
sometimes obscured because the new system is so 
easily solved that one is led directly to the solution 
without noticing the intermediate step." 

The dominant systems can be used for direct com- 
putation of steady states and relaxation dynamics, 
especially when kinetic information is incomplete, 
for design of experiments and mining of experimen- 
tal data, and could serve as a robust first approxi- 
mation in perturbation theory or for precondition- 
ing. They can be used to answer an important ques- 
tion: given a network model, which are its critical 
parameters? Many of the parameters of the initial 
model are no longer present in the dominant sys- 
tem: these parameters are non-critical. Parameters 
of dominant subsystems indicate putative targets to 
change the behavior of the large network. 

Most of reaction networks arc nonlinear, it is 
nevertheless useful to have an efficient algorithm for 
solving linear problems. First, nonlinear systems 
often include linear subsystems, containing reac- 
tions that are (pseudo)monomolecular with respect 
to species internal to the subsystem (at most one 
internal species is reactant and at most one is prod- 
uct). Second, for binary reactions A + B — > if 
concentrations of species A and B (ca,Cb) are well 
separated, say ca 3> cb then we can consider this 
reaction as B — > ... with rate constant proportional 
to ca which is practically constant, because its rel- 
ative changes arc small in comparison to relative 
changes of eg. We can assume that this condition 
is satisfied for all but a small fraction of genuinely 
nonlinear reactions (the set of nonlinear reactions 
changes in time but remains small). Under such 
an assumption, nonlinear behavior can be approxi- 
mated as a sequence of such systems, followed one 
each other in a sequence of "phase transitions" . In 
these transitions, the order relation between some of 
species concentrations changes. Some applications 
of this approach to systems biology are presented by 
Radulescu, Gorban, Zinovyev & Lilienbaum (2008). 
The idea of controllable linearization "by excess" of 
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some reagents is in the background of the efficient 
experimental technique of Temporal Analysis of 
Products (TAP), which allows to decipher detailed 
mechanisms of catalytic reactions (Yablonsky, Olea, 
& Marin (2003)). 

In chemical kinetics various fundamental ideas 
about asymptotical analysis were developed 
(Klonowski (1983)): quasicqiulibrium asymptotic 
(QE), quasi steady-state asymptotic (QSS), lump- 
ing, and the idea of limiting step. 

Most of the works on nonequilibrium thermody- 
namics deal with the QE approximations and correc- 
tions to them, or with applications of these approx- 
imations (with or without corrections). There are 
two basic formulation of the QE approximation: the 
thermodynamic approach, based on entropy max- 
imum, or the kinetic formulation, based on selec- 
tion of fast reversible reactions. The very first use of 
the entropy maximization dates back to the classical 
work of Gibbs (1902), but it was first claimed for a 
principle of informational statistical thermodynam- 
ics by Jaynes (1963). A very general discussion of 
the maximum entropy principle with applications to 
dissipative kinetics is given in the review by Balian, 
Alhassid & Reinhardt (1986). Corrections of QE ap- 
proximation with applications to physical and chem- 
ical kinetics were developed by Gorban, Karlin, Ilg, 
& Ottingcr (2001); Gorban & Karlin (2005). 

QSS was proposed by Bodenstein (1913) and was 
elaborated into an important tool for analysis of 
chemical reaction mechanism and kinetics (Semenov 
(1939); Christiansen (1953); Helfferich (1989)). The 
classical QSS is based on the relative smallness of 
concentrations of some of "active" reagents (radi- 
cals, substrate-enzyme complexes or active compo- 
nents on the catalyst surface) (Aris (1965); Segel & 
Slemrod (1989)). 

Lumping analysis aims to combine reagents into 
"quasicomponcnts" for dimension reduction (Wei & 
Kuo (1969); Kuo & Wei (1969); Li & Rabitz (1989); 
Toth, Li, Rabitz, & Tomlin (1997). 

The concept of limiting step gives the limit simpli- 
fication: the whole network behaves as a single step. 
This is the most popular approach for model simpli- 
fication in chemical kinetics and in many areas be- 
yond kinetics. In the form of a bottleneck approach 
this approximation is very popular from traffic man- 
agement to computer programming and communi- 
cation networks. The proposed asymptotic analysis 
can be considered as a wide extension of the classical 
idea of limiting step (Gorban & Radulescu (2008)). 

The structure of the paper is as follows. In Sec. 2 



we introduce basic notions and notations. We con- 
sider thermodynamic restrictions on the reaction 
rate constants and demonstrate how appear systems 
with arbitrary constants (as subsystems of more de- 
tailed models). For linear networks, the main theo- 
rems which connect ergodic properties with topol- 
ogy of network, are reminded. Four basic ideas of 
model reduction in chemical kinetics are described: 
QE, QSS, lumping analysis and limiting steps. 

In Sec. 3, we introduce the dominant system for a 
simple irreversible catalytic cycle with limiting step. 
This is just a chain of reactions which appears after 
deletion the limiting step from the cycle. Even for 
such simple examples several new observation arc 
presented: 

- The relaxation time for a cycle with limiting step 
is inverse second reaction rate constant; 

- For chains of reactions with well separated rate 
constants left eigenvectors have coordinates close 
to or 1 , and right eigenvectors have coordinates 
close to or ±1. 

For general reaction networks instead of linear 
chains appear general acyclic non-branching net- 
works. For them we also provide explicit formulas 
for eigenvectors and their 0, ±1 asymptotics for 
well-separated constants (Sec. 4). In (Sec. 5) the 
main algorithm is presented. Sec. 6 is devoted to a 
simple demonstration of the algorithm application. 
In Sec. 7, we briefly discuss further corrections to 
dominant systems. The estimates of accuracy are 
given in Appendix. 

2. Main Asymptotic Ideas in Chemical 
Kinetics 

2.1. Chemical Reaction Networks 

To define a chemical reaction network, we have to 
introduce: 

- a list of components (species); 

- a list of elementary reactions; 

- a kinetic law of elementary reactions. 

The list of components is just a list of symbols (la- 
bels) Ai,...A n . Each elementary reaction is repre- 
sented by its stoichiometric equation 

i si 

where s enumerates the elementary reaction, and 
the non-negative integers a s i, /? s j are the stoichio- 
metric coefficients. A stoichiomentric vector 7„ with 
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coordinates -y si — j3 si — a si is associated with each 
elementary reaction. 

For analysis of closed chemical systems with de- 
tailed balance it is usual practice to group reactions 
in pairs, direct and inverse reactions together, but 
in more general settings this is not convenient. 

A non-negative real extensive variable Ni > 0, 
amount of A4 , is associated with each component Ai . 
It measures "the number of particles of that species" 
(in particles, or in moles). The concentration of A, 
is an intensive variable: Cj = Ni/V, where V is vol- 
ume. It is necessary to stress, that in many prac- 
tically important cases the extensive variable V is 
neither constant, nor the same for all components 
Ai. For more details see, for example the book of 
Yablonskii, Bykov, Gorban, & Elokhin (1991). For 
simplicity, we will consider systems with one con- 
stant volume and under constant temperature, but 
it is necessary always keep in mind the possibility 
to return to general equations. For that conditions, 
the kinetic equations have the following form 



dc 



— =J2w s (c,T) Ts + 



(2) 



where v is the vector of external fluxes normalized 
to unit volume. It may be useful to represent exter- 
nal fluxes as elementary reactions by introduction of 
new component together with incoming and out- 
going reactions — > Ai and Ai — > 0. 

The most popular kinetic law of elementary reac- 
tions is the mass action law for perfect systems: 



w, 



(c,T) = k a (T)Y[c?'*, 



(3) 



where "kinetic constant" k s (T) depends on temper- 
ature T. More general kinetic law, which can be used 
for most of non-ideal (non-perfect) systems is 



w s (c, T) 



fs exp a siVi 



(4) 



where R is the universal gas constant, /ij is the chem- 



ical potential, /x, 



dF(N.T.V) _ dG(N,T,P) 



F is the 



dNi dNi 

Helmgoltz free energy, G is the Gibbs energy (free 
enthalpy), P is pressure and (p s > is an intensive 
variable, kinetic factor, which can depend on any set 
of intensive variables, first of all, on T. 

Chemical thermodynamics (Prigogine & Defay 
(1954)) provides tools of choice for stability analy- 
sis of reaction networks (Procaccia & Ross (1977)) 
and chemical reactors (Aris (1965)). The laws of 
thermodynamics have been used for analyzing of 
structural stability of process systems by Hangos, 



Bokor, & Szederkenyi (2004). In general reaction 
network coefficients k s (3) or ip 3 (4) are not inde- 
pendent. In order to respect the second law of ther- 
modynamics, they should satisfy some equations 
and inequalities. The most famous sufficient condi- 
tion gives the principle of detailed balance. Let us 
group the elementary reactions in pairs, direct and 
inverse reactions, and mark the variables for direct 
reactions by superscript +, and for inverse reac- 
tions by — . Then the principle of detailed balance 
for general kinetics (4) reads: 

„+ _ (5) 



(Feinberg (1972)). For the isothermal mass action 
law the principle of detailed balance can be formu- 
lated as follows: there exists a strictly positive point 
c* of detailed balance, at this point 

w t( c *) = w s( c *) (6) 

for all s. This is, essentially, the same principle: if we 
substitute in the general reaction rate (4) the frac- 
tion HijPT by ln(cj/c*), then we will get the mass 
action law, and ipf = ipj. The principle of detailed 
balance is closely related to the microreversibility 
and Onsager relations. 

More general condition was invented by Stueck- 
elberg (1952) for the Boltzmann equation. He pro- 
duced them from the S'-matrix unitarity (the quan- 
tum complete probability formula). For the general 
law (4) without direct-inverse reactions grouping for 
any state the following identity holds: 

Y <Ps exp ^-^ a sitnj 

= Y exp ^-J^ Y PsiVi^J ■ 

Even more general condition which guarantees the 
second law and has clear microscopic sense (the com- 
plete probability does not increase) was obtained by 
Gorban (1984): for any state 



(7) 



s \ i / 

> y <Ps cx p Y P'W^J 



(8) 



To obtain formulas for the isothermal mass action 
law, it is sufficient just to apply the general law (4) 
with constant ip 3 to the perfect free energy F = 
RTj2i c i(^ nc i + Mio) with constant /x,o- More de- 
tailed analysis was presented, by Gorban (1984). 
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In any case, reaction constants are dependent, 
and this dependence guarantees stability of equilib- 
rium and existence of global thermodynamic Lya- 
punov functions for closed systems (2) with v = 
0. Nevertheless, we often study equations for such 
systems with oscillations, bifurcations, chaos, and 
other effects, which are impossible in systems with 
global Lyapunov function. Usually this means that 
we study a subsystem of a large system, where some 
of concentrations do not change because they are 
stabilized by external fluxes or by a large external 
reservoir. These constant (or very slow) concentra- 
tions are included into new reaction constants, and 
after this redefinition they can loose any thermody- 
namic property. 

2.2. Linear Networks and Ergodicity 

In this Sec, we consider a general network of 
linear reactions. This network is represented as a 
directed graph (digraph) (Temkin, Zeigarnik, & 
Bonchev (1996)): vertices correspond to compo- 
nents Ai, edges correspond to reactions Aj — > Aj 
with kinetic constants kji > 0. For each vertex, Ai, 
a positive real variable Cj (concentration) is defined. 
A basis vector e l corresponds to Ai with compo- 
nents e l j = 5ij, where Sij is the Kronecker delta. 
The kinetic equation for the system is 

— j£- = y^(fcjjCj — kjid), (9) 
3 

or in vector form: c = Kc. We don't assume any spe- 
cial relation between constants, and consider them 
as independent quantities. The thermodynamic re- 
strictions on constants are not applicable here be- 
cause, in general, we study pseudomonomolecular 
systems which are subsystems of larger nonlinear 
systems and don't represent by themselves closed 
monomolecular systems. 

For any network of linear reactions the matrix of 
kinetic coefficients K has the following properties: 

- non-diagonal elements of K are non-negative; 

- diagonal elements of K are non-positive; 

- elements in each column of K have zero sum. 
For any K with these properties there exists a net- 
work of linear reactions with kinetic equation c = 
Kc. This family of matrices coincide with the family 
of generators of finite Markov chains, and this class 
of kinetic equations coincide with the class of inverse 
Kolmogorov's equations or master equations for the 



finite Markov chains in continuous time (Meyn & 
Tweedie (2009); Meyn (2007)). 

A linear conservation law is a linear function de- 
fined on the concentrations 6(c) = X^kjCj, whose 
value is preserved by the dynamics (9). The conser- 
vation laws coefficient vectors bi are left eigenvectors 
of the matrix K corresponding to the zero eigen- 
value. The set of all the conservation laws forms the 
left kernel of the matrix K . Equation (9) always has 
a linear conservation law: b°(c) = X)i c * = const. 
If there is no other independent linear conservation 
law, then the system is weakly ergodic. 

A set E is positively invariant with respect to ki- 
netic equations (9), if any solution c(t) that starts 
in E at time to ( c (^o) £ E) belongs to E for t > to 
(c(t) e E if t > t ). It is straightforward to check 
that the standard simplex £ = {c | a > 0, Cj = 
1} is positively invariant set for kinetic equation (9): 
just to check that if Cj = for some i, and all Cj > 
then Cj > 0. This simple fact immediately implies 
the following properties of K: 

- All eigenvalues A of K have non-positive real 
parts, ReX < 0, because solutions cannot leave £ 
in positive time; 

- If ReX = then A = 0, because intersection of 
S with any plane is a polygon, and a polygon 
cannot be invariant with respect to rotations to 
sufficiently small angles; 

- The Jordan cell of K that corresponds to zero 
eigenvalue is diagonal - because all solutions 
should be bounded in £ for positive time. 

- The shift in time operator exp(Kt) is a contrac- 
tion in the l\ norm for t > 0. 

- For weakly ergodic systems there exists such a 
monotonically decreasing function 5(t) (t > 0, < 
5(t) < 1, S(t) —5- when t — > 00) that for any two 
solutions of (9) c(t),c'{t) e S 

^|cj(t)-cU*)|<^)^|cj(0)-cK0)|. (10) 

i i 

The ergodicity coefficient S(t) was introduced by 
Dobrushin (1956) (see also a book by Seneta (1981)). 
It can be estimated using the structure of the net- 
work graph (Gorban, Bykov & Yablonskii (1986); 
Meyn (2007)). 

Two vertices are called adjacent if they share a 
common edge. A path is a sequence of adjacent ver- 
tices. A graph is connected if any two of its vertices 
are linked by a path. A maximal connected sub- 
graph of graph G is called a connected component of 
G. Every graph can be decomposed into connected 
components. 
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A directed path is a sequence of adjacent edges 
where each step goes in direction of an edge. A ver- 
tex A is reachable from a vertex B, if there exists a 
directed path from B to A. 

A nonempty set V of graph vertices forms a sink, if 
there are no directed edges from Ai G V to any Aj £ 
V. For example, in the reaction graph A\ A2 — > 
A3 the one- vertex sets {Ai} and {A3} are sinks. 
A sink is minimal if it does not contain a strictly 
smaller sink. In the previous example, {Ai}, {A3} 
are minimal sinks. Minimal sinks are also called er- 
godic components. 

A digraph is strongly connected, if every vertex A 
is reachable from any other vertex B. Ergodic com- 
ponents are maximal strongly connected subgraphs 
of the graph, but inverse is not true: there may ex- 
ist maximal strongly connected subgraphs that have 
outgoing edges and, therefore, are not sinks. 

The weak ergodicity of the network follows from 
its topological properties. 

Theorem 1. The following properties are equiv- 
alent (and each one of them can be used as an alter- 
native definition of weak ergodicity): 

(i) There exist the only independent linear con- 
servation law for kinetic equations (9) (this is 
b°(c) — J2i c i — const). 

(ii) For any normalized initial state c(0) (b°(c) = 
1) there exists a limit state 

c* = lim vcp(Kt) c(0) 

t— >oo 

that is the same for all normalized initial con- 
ditions: For all c, 

lim exp(Kt)c = b°(c)c*. 

i— >oo 

(hi) For each two vertices Ai, Aj (i ^ j) we can find 
such a vertex A^ that is reachable both from 
Ai and from Aj . This means that the following 
structure exists: 

Ai ->■ ... ->■ A k <- ... <- Aj. 

One of the paths can be degenerated: it may 

be i = k or j = k. 
(iv) The network has only one minimal sink (one 

ergodic component). □ 
The proof of this theorem could be extracted from 
detailed books about Markov chains and networks 
(Meyn (2007); Van Mieghem (2006)). In its present 
form it was published by Gorban, Bykov & Yablon- 
skii (1986) with explicit estimations of ergodicity 
coefficients. 

For every monomolecular kinetic system, the 
maximal number of independent linear conserva- 



tion laws (i.e. the geometric multiplicity of the zero 
eigenvalue of the matrix K) is equal to the maximal 
number of disjoint ergodic components (minimal 
sinks) . 

2.3. Quasi- equilibrium ( QE) or Fast Equilibrium 

Quasi-cquilibrium approximation uses the as- 
sumption that a group of reactions is much faster 
than other and goes fast to its equilibrium. We use 
below superscripts ' f ' and ls ' to distinguish fast and 
slow reactions. A small parameter appears in the 
following form 

A„ 1 

-= £ <(c,T) 7 : + - Y, <(c,T) 7 l, 

(7, slow fast 

To separate variables, we have to study the spaces 
of linear conservation law of the initial system (11) 
and of the fast subsystem 

fast 

If they coincide, then the fast subsystem just dom- 
inates, and there is no fast-slow separation for 
variables (all variables are either fast, or constant). 
But if there exist additional linearly independent 
linear conservation laws for the fast system, then 
let us introduce new variables: linear functions 
b 1 (c) , . . .b n (c) , where b 1 (c) , ...b m (c) is the basis of the 
linear conservation laws for the initial system, and 
fe 1 (c), ...b' m+l (c) is the basis of the linear conservation 
laws for the fast subsystem. Then b m+l+1 (c), ...b n (c) 
are fast variables, b m+1 (c), ...b' m+l (c) are slow vari- 
ables, and & 1 (c), ...b m (c) are constant. The quasi- 
equilibrium manifold is given by the equations 
to* (c, T)"f* — and for small e it serves as an 
approximation to a slow manifold. In the old and 
standard approach it is assumed that system (11) as 
well as system of fast reactions satisfies the thermo- 
dynamic restrictions, and the quasi-equilibrium is 
just a partial thermodynamic equilibrium, and could 
be defined by conditional extremum of thermody- 
namic functions. This guarantees global stability of 
fast subsystems and all the classical singular per- 
turbation theory like Tikhonov theorem could be 
applied. 

Recently, Vora & Daoutidis (2001) took notice 
that this type of reasoning does not require clas- 
sical thermodynamic restrictions on constants. For 
example, let us consider the mass action law ki- 
netics and group the reactions in pairs, direct and 
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inverse reactions. If the set of stoichiometric vec- 
tors for fast reactions is linearly independent, then 
for this system the detailed balance principle holds 
(obviously), and it demonstrates the "thermody- 
namic behaviour" without connection to classical 
thermodynamics. This case of "stoichiometrically 
independent fast reactions" can be generalized for 
irreversible reactions too (Vora& Daoutidis (2001)). 
For such fast system the quasiequilbrium manifold 
has the same nice properties as for thermodynamic 
partial equilibrium, and approximates slow dynam- 
ics for sufficiently small e. 

There are other classes of mass action law sub- 
systems with such a "quasi-thermodynamic" be- 
haviour, which depends on structure, but not on 
constants. For example, any system of reactions 
without interactions has such a property (Gorban, 
Bykov, & Yablonskii (1986)). These reactions have 
the form ctAi — > ^2 any linear reaction are al- 
lowed, as well as reactions like 2Ai — > Aj + 
3Ai — > Aj + Ak + Ai, etc. All such fast subsystems 
can serve for quasi-equilibrium approximation, be- 
cause for them dynamics is globally stable. 

Quasi-equilibrium manifold approximates expo- 
nentially attractive slow manifold and is used in 
many areas of kinetics either as initial approxima- 
tion for slow motion, or just by itself (more discus- 
sion and further references are presented by Gorban 
& Karlin (2005)). 

2.4. Quasi Steady- State ( QSS) or Fast Species 

The quasi steady-state (or pseudo steady state) 
assumption was invented in chemistry for descrip- 
tion of systems with radicals or catalysts. In the 
most usual version the species are split in two groups 
with concentration vectors c s ( "slow" or basic com- 
ponents) and c f ("fast intermediates"). For catalytic 
reactions there is additional balance for c f , amount 
of catalyst, usually it is just a sum 6 f = J2 i c\. The 
amount of the fast intermediates is assumed much 
smaller than the amount of the basic components, 
but the reaction rates are of the same order, or 
even the same (both intermediates and slow compo- 
nents participate in the same reactions). This is the 
source of a small parameter in the system. Let us 
scale the concentrations c f and c s to the compatible 
amounts. After that, the fast and slow time appear 
and we could write c s = W s (c s , c f ), c f = (c s , c f ), 
where e is small parameter, and functions M /s ,VF f 
are bounded and have bounded derivatives (are "of 



the same order" ) . We can apply the standard singu- 
lar perturbation techniques. If dynamics of fast com- 
ponents under given values of slow concentrations 
is stable, then the slow attractive manifold exists, 
and its zero approximation is given by the system 
of equations W f (c s , c f ) = 0. Bifurcations in fast sys- 
tem correspond to critical effects, including ignition 
and explosion. 

This scheme was analyzed many times with plenty 
of details, examples, and some complications. Ex- 
haustive case study of the simplest enzyme reaction 
was provided by Segel & Slemrod (1989) . For het- 
erogenous catalytic reactions, the book by Yablon- 
skii, Bykov, Gorban, & Elokhin (1991) gives analysis 
of scaling of fast intermediates (there are many kinds 
of possible scaling). In the context of the Computa- 
tional Singular Perturbation (CSP) approach, Lam 
(1993) and Lam & Goussis (1994) developed concept 
of the CSP radicals. Gorban & Karlin (2003, 2005) 
considered QSS as initial approximation for slow in- 
variant manifold. Analysis of the error of the QSS 
was provided by Turanyi, Tomlin, & Pilling (1993). 

The QE approximation is also extremely popular 
and useful. It has simpler dynamical properties (re- 
spects thermodynamics, for example, and gives no 
critical effects in fast subsystems of closed systems) . 
Nevertheless, neither radicals in combustion, nor in- 
termediates in catalytic kinetics are, in general, close 
to quasi-equilibrium. They are just present in much 
smaller amount, and when this amount grows, then 
the QSS approximation fails. 

The simplest demonstration of these two approx- 
imation gives the simple reaction: S + E -H- SE — > 
P + E with reaction rate constants kf and k 2 . The 
only possible quasi-equilibrium appears when the 
first equilibrium is fast: kf = K ± /e. The correspond- 
ing slow variable is C s = cs + cse, oe = ce + cse — 
const. For the QE manifold we get a quadratic equa- 
tion jtC SE = c s c E = (C s - c SE ){b E - c S e)- This 
equation gives the explicit dependence cse{C s ), and 
the slow equation reads C s = —k 2 csE{C s ), C s + 
cp = bs = const. 

For the QSS approximation of this reaction ki- 
netics, under assumption bE <C bs, we have fast in- 
termediates E and SE. For the QSS manifold there 
is a linear equation k^csCE — k^csE — k 2 csE = 0, 
which gives us the explicit expression for cse{cs) : 
cse = k^csbE/(kiCs + fcf + k 2 ) (the standard 
Michaelis-Menten formula). The slow kinetics reads 
cs = -kiCs(b E -csE{cs)) + kiC S E(cs)- The differ- 
ence between the QSS and the QE in this example 
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is obvious. 

The terminology is not rigorous, and often QSS is 
used for all singular perturbed systems, and QE is 
applied only for the thermodynamic exclusion of fast 
variables by the maximum entropy (or minimum of 
free energy, or extremum of another relevant ther- 
modynamic function) principle (MaxEnt). This ter- 
minological convention may be convenient. Never- 
theless, without any relation to terminology, the dif- 
ference between these two types of introduction of a 
small parameter is huge. There exists plenty of gen- 
eralizations of these approaches, which aim to con- 
struct a slow and (almost) invariant manifold, and to 
approximate fast motion as well. The following ref- 
erences can give a first impression about these meth- 
ods: Method of Invariant Manifolds (MIM) (Roussel 
& Fraser (1991); Gorban & Karlin (2005), Method of 
Invariant Grids (MIG) , a discrete analogue of invari- 
ant manifolds (Gorban, Karlin, & Zinovyev (2004)), 
Computational Singular Perturbations (CSP) (Lam 
(1993); Lam & Goussis (1994); Zagaris, Kaper, & 
Kaper (2004)) Intrinsic Low-Dimensional Manifolds 
(ILDM) by Maas, & Pope (1992), developed further 
in series of works by Bykov, Goldfarb, Gol'dshtein, & 
Maas, U. (2006)), methods based on the Lyapunov 
auxiliary theorem (Kazantzis & Kravaris (2006)). 

2.5. Lumping Analysis 

Wei & Prater (1962) demonstrated that for 
(pseudo)monomolecular systems there exist linear 
combinations of concentrations which evolve in time 
independently. These linear combinations (quasi- 
components) correspond to the left eigenvectors of 
kinetic matrix: if IK = XI then d(l,c)/dt = (l,c)X, 
where the standard inner product (I, c) is concentra- 
tion of a quasicomponent. They also demonstrated 
how to find these quasicomponents in a properly 
organized experiment. 

This observation gave rise to a question: how 
to lump components into proper quasicomponents 
to guarantee the autonomous dynamics of the 
quasicomponents with appropriate accuracy. Wei 
and Kuo studied conditions for exact (Wei & Kuo 
(1969)) and approximate (Kuo & Wei (1969)) lump- 
ing in monomolecular and pseudomonomolecular 
systems. They demonstrated that under certain 
conditions large monomolecular system could be 
well-modelled by lower-order system. 

More recently, sensitivity analysis and Lie group 
approach were applied to lumping analysis (Li & 



Rabitz (1989); Toth, Li, Rabitz, & Tomlin (1997)), 
and more general nonlinear forms of lumped con- 
centrations are used (for example, concentration of 
quasicomponents could be rational function of c) . 

Hutchinson & Luss (1970) studied lumping- 
analysis of mixtures with many parallel first order 
reactions. Farkas (1999) generalized these results 
and characterized those lumping schemes which 
preserve the kinetic structure of the original system. 
Coxson & Bischoff (1987) placed lumping analysis 
in the linear systems theory and demonstrated the 
relationships between lumpability and the concepts 
of observability, controllability and minimal real- 
ization. Djouad & Sportisse (2002) considered the 
lumping procedures as efficient techniques leading 
to nonstiff systems and demonstrated efficiency 
of developed algorithm on kinetic models of at- 
mospheric chemistry. Lin, Leibovici & Jorgensen 
(2008) formulated an optimal lumping problem as 
a mixed integer nonlinear programming (MINLP) 
and demonstrated that it can be efficiently solved 
with a stochastic optimization method, Tabu Search 
(TS) algorithm. 

The power of lumping using a time-scale based 
approach was demonstrated by Whitchousc, Tom- 
lin, & Pilling (2004). This computationally cheap 
approach combines ideas of sensitivity analysis with 
simple and useful grouping of species with similar 
lifetimes and similar topological properties caused 
by connections of the species in the reaction net- 
works. The lumped concentrations in this approach 
are simply sums of concentrations in groups. For ex- 
ample, species with similar composition and func- 
tionalities could be lumped into one single represen- 
tative species (Pepiot-Desjardins & Pitsch (2008)). 

Lumping analysis based both on mathematical 
arguments and fundamental physical and chemi- 
cal properties of the components is now one of the 
main tools for model reduction in highly multicom- 
ponent systems, such as the hydrocarbon mixture 
in petroleum chemistry (Zavala & Rodriguez & 
Vargas- Villamil (2004)) or biochemical networks 
in systems biology (Maria (2006)). The optimal 
solution of lumping problem often requires the ex- 
haustive search, and instead of them various heuris- 
tics are used to avoid combinatorial explosion. For 
the lumping analysis of the systems biology mod- 
els Dokoumetzidis & Aarons (2009) developed a 
heuristic greedy search strategy which allowed them 
to avoid the exhaustive search of proper lumped 
components. 

Procedures of lumping analysis form a part of gen- 



8 



eral algebra of model building and model simplifi- 
cation transformations. Hangos & Cameron (2001) 
applied formal methods of computer science and ar- 
tificial intelligence for analysis of this algebra. In 
particular, a formal method for defining syntax and 
semantics of process models has been proposed. 

The modern systems and control theory provides 
efficient tools for lumping-analysis. The so-called 
balanced model reduction was invented in late 1970s 
(Moore (1981)). For a linear system a set of "target 
variables" is selected. The dimension of the system n 
is large, while the number of the target variables, for 
example, inputs m and outputs p, usually satisfies 
m,p <C n. The balanced model reduction problem can 
be stated as follows (Gugercin & Antoulas (2004)): 
find a reduced order system such that the following 
properties are satisfied: 

(i) The approximation error in the target vari- 
ables is small, and there exists a global error 
bound. 

(ii) System properties, like stability and passivity, 
are preserved. 

(iii) The procedure is computationally efficient. 
In large dimensions, special efforts are needed to re- 
solve the accuracy /efficiency dilemma and to find 
efficiently the approximate solution of the model re- 
duction problem (Antoulas & Sorensen (2002)). 

Various methods for balanced truncation are de- 
veloped: Lyapunov balancing, stochastic balancing, 
bounded real balancing, positive real balancing, and 
frequency weighted balancing (Gugercin & Antoulas 
(2004)). Nonlinear generalizations are proposed as 
well (Lall, Marsden & Glavaki (2002); Condon & 
Ivanov (2004)). 

2.6. Limiting Steps 

In the IUPAC Compendium of Chemical Ter- 
minology (2007) one can find a definition of lim- 
iting steps. Rate-controlling step (2007): "A rate- 
controlling (rate-determining or rate-limiting) step 
in a reaction occurring by a composite reaction se- 
quence is an elementary reaction the rate constant 
for which exerts a strong effect - stronger than that 
of any other rate constant - on the overall rate." 

Let us complement this definition by additional 
comment: usually when people are talking about 
limiting step they expect significantly more: there 
exists a rate constant which exerts such a strong ef- 
fect on the overall rate that the effect of all other rate 
constants together is significantly smaller. For the 



IUPAC Compendium definition a rate-controlling 
step always exists, because among the control func- 
tions generically exists the biggest one. On the con- 
trary, for the notion of limiting step that is used 
in practice, there exists a difference between sys- 
tems with limiting step and systems without limit- 
ing step. 

During XX century, the concept of the limiting 
step was revised several times. First simple idea of a 
"narrow place" (the least conductive step) could be 
applied without adaptation only to a simple cycle 
or a chain of irreversible steps that are of the first 
order (see Chap. 16 of the book Johnston (1966) or 
the paper by Boyd (1978)). When researchers try to 
apply this idea in more general situations they meet 
various difficulties such as: 

- Some reactions have to be "pseudomonomolecu- 
lar." Their constants depend on concentrations 
of outer components, and are constant only un- 
der condition that these outer components are 
present in constant concentrations, or change suf- 
ficiently slow (i.e. are present in significantly big- 
ger amount). 

- Even under fixed or slow outer components con- 
centration, the simple "narrow place" behaviour 
could be spoiled by branching or by reverse reac- 
tions. The simplest example is given by the cycle: 
A\ o A2 — > A3 — > A\. Even if the constant of 
the last step A 3 — > A\ is the smallest one, the 
stationary rate may be much smaller than k^b 
(where b is the overall balance of concentrations, 
b = c\ + c 2 + c 3 ), if the constant of the reverse 
reaction A 2 — > Ai is sufficiently big. 

In a series of papers, Northrop (1981, 2001) clearly 
explained these difficulties and suggested that the 
concept of rate-limiting step is "outmoded" . Nev- 
ertheless, the main idea of limiting is so attractive 
that Northrop's arguments stimulated the search for 
modification and improvement of the main concept. 

Ray (1983) proposed the use of sensitivity analy- 
sis. He considered cycles of reversible reactions and 
suggested a definition: The rate-limiting step in a 
reaction sequence is that forward step for which a 
change of its rate constant produces the largest effect 
on the overall rate. 

Ray's approach was revised by Brown & Cooper 
(1993) from the system control analysis point of 
view (see the book of Cornish-Bowden & Cardenas 
(1990)). They stress again that there is no unique 
rate-limiting step specific for an enzyme, and this 
step, even if it exists, depends on substrate, product 
and effector concentrations. 
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Near critical conditions the critical simplification 
appears, which is also a type of limitation, because 
some reactions become critically important (Yablon- 
sky, Mareels, & Lazman (2003)) 

Two classical examples of limiting steps demon- 
strate us the chain of linear reaction and the linear 
catalytic cycle, when they include a reaction which 
is significantly slower, than other reactions. 

A linear chain of reactions, A\ — > A 2 — > ■■■A n , 
with reaction rate constants fcj (for A, — > Aj+i), 
gives the first example of limiting steps. Let the 
reaction rate constant k q be the smallest one. Then 
we expect the following behaviour of the reaction 
chain in time scale > l/k q : all the components 
A\, ...A q -i transform fast into A q , and all the com- 
ponents A q+ i, ...A n -i transform fast into A n , only 
two components, A q and A n are present (concentra- 
tions of other components are small) , and the whole 
dynamics in this time scale can be represented by 
a single reaction A q —> A n with reaction rate con- 
stant k q . This picture becomes more exact when k q 
becomes smaller with respect to other constants. 

The catalytic cycle is one of the most important 
substructures that we study in reaction networks. In 
the reduced form the catalytic cycle is a set of linear 
reactions: 



Ai ->■ A 2 ->■ . . . A n 



Al 



Reduced form means that in reality some of these 
reaction are not monomolecular and include some 
other components (not from the list A\, . . . A n ). But 
in the study of the isolated cycle dynamics, concen- 
trations of these components are taken as constant 
and are included into kinetic constants of the cycle 
linear reactions. 

For the constant of elementary reaction Ai — > we 
use the simplified notation fc, because the product 
of this elementary reaction is known, it is for 
i < n and A\ for i = n. The elementary reaction 
rate is u>i = kid, where Cj is the concentration of 
Ai. The kinetic equation is: 



Cj ki—\Ci—\ kiCi 



(12) 



where by definition co — c n , k — k n , and wq — 
w n . In the stationary state (cj = 0), all the Wi are 
equal: Wi — w. This common rate w we call the cycle 
stationary rate, and 

b w 



w 



' 1 c i 



■ ■ ■ k„ 



k'i 



(13) 



where b = J2 i Ci is the conserved quantity for reac- 
tions in constant volume. Let one of the constants, 



&min, be much smaller than others (let it be /c m ; n = 
h > k 

min 

if i^n . (14) 
In this case, in linear approximation w = k n b, 

c n = b 1 1 — ) , and a = b-^- for i 7^ n . 

(15) 

The simplest zero order approximation for the 
steady state gives 



c n = b, Ci = (i ^ n). 



(16) 



This is trivial: all the concentration is collected at 
the starting point of the "narrow place," but may be 
useful as an origin point for various approximation 
procedures. 

So, the stationary rate of a cycle is determined 
by the smallest constant, fc m ; n , if it is much smaller 
than the constants of all other reactions (14): 



w w k min b. 



(17) 



In that case we say that the cycle has a limiting step 
with constant fc m ; n . 

3. Dynamics of Catalytic Cycle with 
Limiting Step 

3.1. Eigenvalues 

There is significant difference between the exam- 
ples of limiting steps for the chain of reactions and 
for irreversible cycle. For the chain, the steady state 
does not depend on nonzero rate constants. It is just 
c n = b, c\ = c 2 = ... = c„_i = 0. The smallest rate 
constant k q gives the smallest positive eigenvalue, 
the relaxation time is r = l/k q . The corresponding 
approximation of eigenmode (right eigenvector) r 1 
has coordinates: r\ = ... = r q _ 1 = 0, r q = 1, r q+1 = 
... = r n _i = 0,r n = —1. This exactly corresponds to 
the statement that the whole dynamics in the time 
scale > 1/kg can be represented by a single reaction 
A q — > A n with reaction rate constant k q . The left 
eigenvector for eigenvalue k q has approximation I 1 
with coordinates l\ = l\ = ... = l q = 1, l q+1 = ... = 
l n = 0. This vector provides the almost exact lump- 
ing on time scale > l/k q . Let us introduce a new 
variable ci ump = J2i l *Ci, i.e. ci ump = ci+c 2 + ... + c q . 
For the time scale > l/k q we can write q ump + c n ~ 
b, dci um p/dt ~ k q c\ um p, dc n /dt ~ ^Qump- 

In the example of a cycle, we approximate the 
steady state, that is, the right eigenvector r° for 
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zero eigenvalue (the left eigenvector is known and 
corresponds to the main linear balance b: I® = 1). In 
the zero-order approximation, this eigenvector has 
coordinates r\ = ... = r°_ x = 0, r" = 1. 

If k n /ki is small for all i < n, then the kinetic 
behaviour of the cycle is determined by a linear chain 
ofn— 1 reactions A\ — > A2 — > ...A n , which we obtain 
after cutting the limiting step. The characteristic 
equation for an irreversible cycle, n™=i(^ + k i) — 
n"=i hi = 0; tends to the characteristic equation for 
the linear chain, A J\7=i + ki ) = ^> wnen — ^ 0. 

The characteristic equation for a cycle with limit- 
ing step {kn/ki <C 1) has one simple zero eigenvalue 
that corresponds to the conservation law Cj = b 
and n — 1 nonzero eigenvalues 

A 4 = -fcj + (i < n). (18) 

where Si -> when X) i<n I s " ^ °- 

A cycle with limiting step (12) has real eigenspec- 
trum and demonstrates monotonic relaxation with- 
out damped oscillations. Of course, without limita- 
tion such oscillations could exist, for example, when 
all ki = k > 0, (i = 1, ...n). 

The relaxation time of a stable linear system (12) 
is, by definition, r = 1/ min{i?e(— A»)} (A 7^ 0). For 
small fc„, r w l/fc T , k T = minjfcj}, (« = 1, ...n — 1). 
In other words, for a cycle with limiting step, k T is 
the second slowest rate constant: fc m j n <C fc T < •••• 

3.2. Eigenvectors for Reaction Chain and for 
Catalytic Cycle with Limiting Step 

Let the irreversible cycle include a limiting step: 
k n <C fcj (i — 1, ...,n — 1) and, in addition, k n <C 
|fcj — fcj I — 1, n — 1, i ^ j), then the eigenvec- 
tors of the kinetic matrix almost coincide with the 
eigenvectors for the linear chain of reactions A\ — > 
A2 — > ...A n , with reaction rate constants fc, (for 
Ai -> Ai+i) (Gorban & Radulescu (2008)). 

The kinetic equation for the linear chain is 



c'i — ki—\Ci — \ kiCi 



(19) 



The coefficient matrix K of these equations is very 
simple. It has nonzero elements only on the main 
diagonal, and one position below. The eigenvalues 
of K are — ki (i — 1, ...n — 1) and 0. The left and 
right eigenvectors for eigenvalue, 1° and r°, are: 

/° = (l,l,...l), r° = (0,0,...0,l), (20) 

all coordinates of 1° are equal to 1, the only nonzero 
coordinate of r° is r n and we represent vector- 
column r° in row. 



Below we use explicit form of K left and right 
eigenvectors. Let vector-column r l and vector-row 
l l be right and left eigenvectors of K for eigenvalue 
— ki. For coordinates of these eigenvectors we use 
notation r] and 1*. Let us choose a normalization 
condition r\ = l\ = 1. It is straightforward to check 
that r) = (j < i) and I) = (j > i), r) +1 = 

k j r j/( k 3 + l - k i) U > i) and l j-l = k j-l l j/( k 3-l - 

kj) (j < i), and 



in j 

i = TT K-i+3-1 . ,i 
i+rn 1 1 , _ , > h- 



n 



ki—j k{ 



■ (21) 



It is convenient to introduce formally fco = 0. Under 
selected normalization condition, the inner product 
of eigenvectors is: ZV- 7 = Sij, where Sij is the Kro- 
necker delta. 

If the rate constants any two constants, ki, kj are 
connected by relation ki 3> kj or ki <C kj (i.e. they 
are well separated), then 



if ki < ki 



l , ■ (22) 

k{— j ki y 0, if k^ 



Hence, |Z*_ m 
that also |r 



1 or \l\_ m \ ~ 0. To demonstrate 

i+m\ ~ 1 or KVml ~ °. wc snift nomina- 
tors in the product (21) on such a way: 



i+m 



k%-\-r 



ki. 



n 



ki 



here is 



Exactly as in (22), each multiplier 
either almost 1 or almost 0, and k . + ki _ k . is either 
almost or almost — 1. In this zero-one asymptotics 



I) =1, D 



1 



if fci-j > fcj for all j = 1, ... to, else 



0; 



-1 



if fci+j > fcj for all j = 1, ... m — 1 
and fc l+m < ki, else r| +m « 0. 

(23) 

In this asymptotic (Fig. 1), only two coordinates of 
right eigenvector r % can have nonzero values, r\ = l 
and r\ +m w — 1 where to is the first such positive 
integer that i + m < n and k i+m < ki. Such to 
always exists because k n — 0. For left eigenvector 
l % , l\ w . . . w 1 and Z|_ • w where j > and 
g is the first such positive integer that i — q — I > 
and ki- q -i < ki. It is possible that such g does 
not exist. In that case, all Z| • w 1 for j > 0. It 
is straightforward to check that in this asymptotic 
= S^. 
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i 

Fig. 1. Graphical representation of eigenvectors approxima- 
tion for the linear chain of reactions with well separated 
constants. To find the left (Z) and right (r) eigenvectors for 
eigenvalue k it is necessary to delete from the chain all the 
reactions with the rate constants < k (dashed lines) and to 
find the maximal connected interval, where the reaction with 
constant k (bold arrow) is situated. The right eigenvector r 
has coordinate 1 for the vertex, which is the beginning of the 
reaction with constant k, and coordinate —1 for the vertex, 
which is end of the interval in the direction of reactions. The 
left eigenvector I has coordinate 1 for the beginning of the 
reaction with constant k and for all preceding vertices from 
the connected interval. All other coordinates of r and I are 
zero. 

The simplest example gives the order k\ ^> k 2 ~> 
... » fc„_i: « 1 for j > 0, r\ = 1, rj +1 » -1 
and all other coordinates of eigenvectors are close to 
zero. For the inverse order, k\ <C k 2 <C ... <C fc„-i, 
l\ = 1, r\ = 1, r l n « — 1 and all other coordinates of 
eigenvectors are close to zero. 

For less trivial example, let us find the asymptotic 
of left and right eigenvectors for a chain of reactions: 

A 1 hA 2 hA i hA i hA 5 hA 6 , 

where the upper index marks the order of rate con- 
stants: ii4 » ^5 » ^ » ^ » k\ (ki is the rate 
constant of reaction Ai — > ...). 

For left eigenvectors, rows l l , we have the following 
asymptotics: 

I 1 « (1,0,0,0,0,0), Z 2 « (0,1,0,0,0,0), 

Z 3 «(0,1,1,0,0,0),Z 4 «(0,0,0,1,0,0), (24) 

Z 5 « (0,0,0,1,1,0). 

For right eigenvectors, columns r l , we have the 
following asymptotics (we write vector-columns in 
rows): 

r 1 ^ (1,0,0,0,0,-1), r 2 « (0,1,-1,0,0,0), 

r 3 w (0, 0, 1,0, 0, -1), r 4 « (0, 0, 0, 1, -1, 0), (25) 

r 5 « (0,0,0,0,1,-1). 

The corresponding approximation to the general so- 
lution of the kinetic equations is: 

n-l 

c(t) = (1°, c(0))r° + ^(fc(0))r 2 cxp(-M), (26) 

where c(0) is the initial concentration vector, and 
for left and right eigenvectors l % and r l we use their 



zero-one asymptotic. In other words, approximation 
of the left eigenvectors provides us with almost exact 
lumping (for analysis of exact lumping see the paper 
by Li & Rabitz (1989)) . 

4. Acyclic Non-branching Network: Explicit 
Formulas for Eigenvectors 

So, to analyze asymptotic of eigenvalues and 
eigenvectors for a irreversible cycle, we cut the reac- 
tion with the smallest constant, get a linear chain, 
and analyze the eigenvalues and eigenvectors for 
this chain. For a general multiscale reaction net- 
work (instead of a cycle) we will come, after some 
surgery, to acyclic non-branching reaction networks 
(instead of a linear chain) . 

For any network without branching, we can sim- 
plify the notation for the kinetic constants, by intro- 
ducing Ki = kji for the only reaction Ai — > Aj, or 
Ki = 0, if there is no such a reaction. Also it is useful 
to introduce a map <j> on the set of vertices: <f>(i) = j, 
if there exist reaction Ai — > Aj , and <j>(i) = i if there 
are no outgoing reactions from the Ai -4- Aj. For 
iterations of the map <fi we use notation <\fl . 

For an acyclic non-branching reaction network, 
for any vertex A; t there is an eigenvalue — Ki and 
the corresponding eigenvector. If Ai is a sink vertex, 
then this eigenvalue is zero. For left and right eigen- 
vectors of K that correspond to Ai we use notations 
l l (vector- row) and r % (vector-column) , correspond- 
ingly. 

Let us suppose that Af is a sink vertex of the 
network. Its associated right and left eigenvectors 
corresponding to the zero eigenvalue are given by: 
Tj = 6ij\ l l j = 1 if and only if (j) q {j) = i for some 
q > 0. 

For nonzero eigenvalues, right eigenvectors will be 
constructed by recurrence starting from the vertex 
Ai and moving in the direction of the flow. The con- 
struction is in opposite direction for left eigenvec- 
tors. 

For right eigenvector r l only coordinates r^ fc ^ 
(k = 0, 1, . . .Ti) could have nonzero values, and 

k 

i _ K fc W i _ TT K <t> 3 W 

- iu ^ AJL _ Kt 

fe-i 

_ Ki TT + 

K^k+l^i) - K l - K l 

(27) 
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Fig. 2. Graphical representation of eigenvectors approxima- 
tion for the acyclic non-branching reaction network with well 
separated constants (compare to Fig. 1). The eigenvalue — k 
corresponds to the reaction Aj — > A^(j\ (bold arrow). To the 
right from Aj are vertices A^qi^ and to the left are those 
Aj, for which there exists such q that 4> q (j) = i. The reac- 
tions with the rate constants < k (dashed lines) are deleted 
from the network. The right and left eigenvectors could have 
nonzero coordinates only for vertices from the maximal con- 
nected subgraph of the presented graph, where the Aj is 
situated. The right eigenvector r has coordinate 1 for Aj 
(beginning of the bold arrow), and coordinate —1 for the 
vertex, which is the minimal in that connected subgraph. 
The left eigenvector I has coordinate 1 for the beginning of 
the reaction with constant k and for all preceding vertices 
from the subgraph. All other coordinates of r and I are zero. 

For left eigenvector l l coordinate lj could have 
nonzero value only if there exists such q > that 
^ 9 (j) = i (this q is unique because the system is 
acyclic): 

Kj Ki K t 

For well separated constants, we can write the 
asymptotic representation explicitly, analogously to 
(23) (Fig. 2). For left eigenvectors, l\ = 1 and I) = 
1 (for i ^ j) if there exists such q that <fi q (j) = i, 
and K^dfj) > ^ for all d — 0, . . . q — 1, else i*- = 0. 
For right eigenvectors, r\ = 1 and r l ^ k ^ = — 1 if 
K^*(i) < Ki and for all positive m < k inequality 
> ^ holds, i.e. k is first such positive inte- 
ger that < Ki (for fixed point A p we use k p = 
0). Vector r l has not more than two nonzero coor- 
dinates. It is straightforward to check that in this 
asymptotic i'r* = <%. 

For example, let us find that asymptotic for a 
branched acyclic system of reactions: 

A 1 ^A 2 hA 3 hA 4 hA^A s , A 6 hA 7 hA 4 

where the upper index marks the order of rate con- 
stants: Kg > K4 > K7 > K5 > K2 > «3 > K\ {Ki is 

the rate constant of reaction Aj, — > ...). 

For zero eigenvalue, the left and right eigenvectors 
are 

Z 8 = (1,1, 1,1, 1,1, 1,1,1), r 8 = (0,0,0,0,0,0,0,1). 



For left eigenvectors, rows l l , that correspond to 
nonzero eigenvalues we have the following asymp- 
totics: 

I 1 w (1,0,0,0,0,0,0,0), I 2 w (0,1,0,0,0,0,0,0), 
I 3 w (0, 1,1, 0,0,0, 0,0), Z 4 « (0,0,0,1,0,0,0,0), 
Z 5 w (0, 0, 0, 1, 1, 1, 1, 0), f « (0, 0, 0, 0, 0, 1, 0, 0). 
Z 7 w (0,0,0,0,0,1,1,0) 

(29) 

For the corresponding right eigenvectors, columns 
r l , we have the following asymptotics (we write 
vector-columns in rows): 

r 1 ^ (1,0, 0,0, 0,0, 0,-1), r 2 w(0, 1,-1, 0,0, 0,0,0), 
r 3 «(0, 0,1,0,0,0,0,-1), r 4 «(0, 0,0, 1,-1, 0,0,0), 
r 5 « (0,0,0, 0,1,0,0,-1), r 6 «(0, 0,0, 0,0, 1,-1,0), 
r 7 «(0, 0,0, 0,-1, 0,1,0). 

(30) 

5. Calculating the Dominant System for a 
Linear Multiscale Network 

5.1. Problem Statement 

We study asymptotical behavior of the transfor- 
mation of the kinetic matrix K to the normal form 
along the lines lnfcjj = when £ — > oo. For al- 
most all direction vectors (%) (outside several hy- 
pcrplanes) there exists a minimal reaction network 
which reaction rate constants are monomials of kij 

(Ylij Hj' 1 where are not obligatory positive num- 
bers) and eigenvectors and eigenvalues approximate 
the eigenvectors and eigenvalues when £ — > oo with 
arbitrary high relative accuracy. We call this mini- 
mal system the dominant system. Existence of dom- 
inant systems is proven by direct construction (this 
Sec.) and estimates of accuracy of approximations 
(Appendix) . 

The dominant systems coincide for vectors (%) 
from some polyhedral cones. Therefore, we don't 
need to study a given value of (9ij) but rather have 
to build these cones together with the correspon- 
dent dominant systems. The following formal rule 
( "assumption of well separated constants" ) allows 
us to simplify this task: if in construction of dom- 
inant systems we need to compare two monomials, 

Mf = Ua Hi and M 9 = Uij fcy* then ^ can al- 
ways state that either Mf 3> M g or Mf <C M g and 
consider the logarithmic hyperplane Mf = M g as a 
boundary between different cones. At the end, we 



13 



Aj ^kji=vnzK l {k li } 
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Fig. 3. Construction of the auxiliary reaction network by 
pruning. For every vertex, it is necessary to leave the out- 
going reaction with maximal reaction rate constant. Other 
reactions should be deleted. 

can join all cones with the same dominant system. 
We are interested in robust asymptotic and do not 
analyze directions ) which belong to the bound- 
ary hyperplanes. This robust asymptotic with well 
separated constants and acyclic dominant systems 
is typical because the exclusive direction vectors be- 
lon to a finite number of hyperplanes. 

There may be other approaches based on (i) the 
Maslov dequantization and idempotent algebras 
(Litvinov & Maslov (2005)), (ii) the limit of log- 
uniform distributions in wide boxes of constants 
under some conditions (Feng, Hooshangi, Chen, 
Li, Weiss, & Rabitz (2004); Gorban & Radulescu 
(2008)), or (iii) on consideration of all possible or- 
derings of all monomials with integer exponents and 
construction of correspondent dominant systems 
(Robbiano (1985) proved that there exists only a fi- 
nal number of such orderings and enumerated all of 
them, see also the book by Greuel & Pfister (2002)). 
They give the same final result but with different 
intermediate steps. 

5.2. Auxiliary Operations 

5.2.1. From Reaction Network to Auxiliary 
Dynamical System 

Let us consider a reaction network W with a given 
structure and fixed ordering of constants. The set of 
vertices of W is A and the set of elementary reactions 
is 1Z. Each reaction from 7Z has the form Ai — > Aj , 
Ai,Aj € A. The corresponding constant is kji. For 
each Ai G A we define Hi = maxjjfcji} and (f>(i) = 
arg m&Xj{kji}. In addition, <p(i) = i if kji — for all 
3- 

The auxiliary discrete dynamical system for the re- 
action network W is the dynamical system $ = $w 
defined by the map <j> on the finite set A. The auxil- 
iary reaction network (Fig. 3) V = Vyv has the same 
set of vertices A and the set of reactions Ai — > A^u-s 
with reaction constants Ki. Auxiliary kinetics is de- 




* * » 




Fig. 4. Decomposition of a discrete dynamical system. 



scribed by c = Kc, where = —K,j^>ij + Kjdi^jy 



5.2.2. Decomposition of Discrete Dynamical 
Systems on Finite Sets 

Discrete dynamical system on a finite set V = 
{Ai,A 2 , . . . A n } is a semigroup 1,0, 2 , where <j) 
is a map <f> : V — > V. Ai € V is a periodic point, 
if (j) (Ai) = Ai for some I > 0; else Ai is a tran- 
sient point. A cycle of period I is a sequence of I 
distinct periodic points A, (j>(A), 4> 2 {A), . . . (j) l ^ 1 (A) 
with <f) l (A) = A. A cycle of period one consists of 
one fixed point, 4>(A) = A. Two cycles, C, C cither 
coincide or have empty intersection. 

The set of periodic points, V p , is always 
nonempty. It is a union of cycles: V p — UjCj. For 
each point A e V there exist such a positive integer 
t{A) and a cycle C(A) = C 3 that <j>i(A) e Cj for 
q > t(A) . In that case we say that A belongs to basin 
of attraction of cycle Cj and use notation Att(Cj) = 
{A | C(A) = C-j}. Of course, Cj C Att{Cj). For dif- 
ferent cycles, Att(Cj)nAtt(Ci) = 0. If A is periodic 
point then t(A) = 0. For transient points t(A) > 0. 

So, the phase space V is divided onto subsets 
Att(Cj) (Fig. 4). Each of these subsets includes 
one cycle (or a fixed point, that is a cycle of length 
1). Sets Att(Cj) are ^-invariant: 4>(Att(Cj)) C 
Att(Cj). The set Att(Cj) \ Cj consist of transient 
points and there exists such positive integer r that 
<ffl(Att(Cj)) = Cj if q > t. 

Discrete dynamical systems on a finite sets corre- 
spond to graphs without branching points. Notice 
that for the graph that represents a discrete dy- 
namic system, attractors are ergodic components, 
while basins are connected components. 



14 



c ^ 

/ * r 
\ / X,p I 

Fi{|. 5. Gluing a cycle with rate constants rcnormalization. 
c^ S are the quasistationary concentrations on the cycle. Af- 
ter gluing, we have to leave the outgoing from A 1 reaction 
with the maximal renormalized rate constant, and delete 
others. 

5.3. Algorithm for Calculating the Dominant 
System 

For this general case, the algorithm consists of 
two main procedures: (i) cycles gluing and (ii) cycles 
restoration and cutting. 

5.3.1. Cycles Gluing 

Let us start from a reaction network W with a 
given structure and fixed ordering of constants. The 
set of vertices of W is A and the set of elementary 
reactions is 1Z. 

If all attractors of the auxiliary dynamic system 
$w are fixed points Af\,Af 2 , ... 6 A, then the aux- 
iliary reaction network is acyclic, and the auxiliary 
kinetics approximates relaxation of the whole net- 
work W. 

In general case, let the system $>v have sev- 
eral attractors that are not fixed points, but cycles 
Ci,C2,... with periods ti,t 2 ,... > 1. By gluing 
these cycles in points, we transform the reaction 
network W into W 1 . The dynamical system $yv 
is transformed into $ 1 . For these new system and 
network, the connection = $ w i persists: $ 1 is 
the auxiliary discrete dj d system for W . 

For each cycle, d, we introduce a new vertex A 1 . 
The new set of vertices, A 1 = A U {A 1 , A 2 , ...} \ 
(UjCj) (we delete cycles Ci and add vertices A 1 ). 

All the reaction A — > B from the initial set 1Z, 
(A, B G A) can be separated into 5 groups: 

(i) both A,B £ \J t C l ; 

(ii) A £ [Jid, but B e C t ; 
(in) A £ d, but B £ U i C l ; 

(iv) AeC„Be Cj, i^i; 

(v) i,BeC, 

Reactions from the first group do not change. Reac- 
tion from the second group transforms into A — > A 1 



(to the whole glued cycle) with the same constant. 
Reaction of the third type changes into A 1 — > B with 
the rate constant renormalization: let the cycle C % 
be the following sequence of reactions A\ — > A 2 — !> 
...A Ti — > A\ , and the reaction rate constant for A; — > 
Ai + i is fcj (k Ti for A Ti — > A\). For the limiting reac- 
tion of the cycle Cj we use notation fcu m j. If A = Aj 
and k is the rate reaction for A — > S, then the new 
reaction A 1 ^ B has the rate constant fcfcii m i/fcj. 
This corresponds to a quasistationary distribution 
on the cycle (15). The new rate constant is smaller 
than the initial one: kku m i/kj < k, because fcn m , < 
kj due to definition of limiting constant. The same 
constant renormalization is necessary for reactions 
of the fourth type. These reactions transform into 
A 1 — > A J . Finally, reactions of the fifth type vanish. 

After we glue all the cycles (Fig. 5) of auxiliary 
dynamical system in the reaction network W, we get 
W 1 . Let us assign W := W 1 , A := A 1 and iterate 
until we obtain an acyclic network and exit. This 
acyclic network is a "forest" and consists of trees 
oriented from leafs to a root. The number of such 
trees coincide with the number of fixed points in the 
final network. 

After gluing we can identify the reactions, which 
will be included into the dominant system. Their 
constants are the critical parameters of the networks. 
The list of these parameters, consists of all reac- 
tion rates of the final acyclic auxiliary network, and 
of the rate constants of the glued cycles, but with- 
out their limiting steps. Some of these parameters 
are rate constants of the initial network, other have 
the monomial structure. Other constants and corre- 
sponding reactions do not participate in the follow- 
ing operations. To form the structure of the domi- 
nant network, we need one more procedure. 

5.3.2. Cycles Restoration and Cutting 

We start the reverse process from the glued net- 
work V m on A m . On a step back, from the set A m to 
A" 1 ^ 1 and so on, some of glued cycles should be re- 
stored and cut. On the qth step we build an acyclic 
reaction network on A m ~ q , the final network is de- 
fined on the initial vertex set and approximates re- 
laxation of W. 

To make one step back from V m let us select the 
vertices of A m that are glued cycles from V™ -1 . Let 
these vertices be A™, A™, .... Each A™ corresponds 
to a glued cycle from V m ~\ A™' 1 -> A™' 1 -> 
...A™- 1 ->• A™' 1 , of the length Ti . We assume that 
the limiting steps in these cycles are A™.' 1 — > A™ -1 . 



15 



h .... 4 

t T 1 ^=> ^ 

\ / ^ 

\ 1 kk^ m l k t 

-. 

' 4 

Fig. 6. The main operation of the cycle surgery: on a step 
back we get a cycle A\ —»...—> A T — ► A\ with the limiting 
step At — ¥ A\ and one outgoing reaction Ai — > Aj. We 
should delete the limiting step, reattach ("recharge") the 
outgoing reaction Aj — s- Aj from Ai to A r and change its 
rate constant k to the rate constant kku m /ki. The new value 
of reaction rate constant is always smaller than the initial 
one: fcfcii m /fci < k if fcij m ^ ki. For this operation only one 
condition k <C fci is necessary (A; should be small with respect 
to reaction Ai — > rate constant, and can exceed any 

other reaction rate constant). 

Let us substitute each vertex A™ in V m by n ver- 
tices A™-\A%-\ ...A™' 1 and add to V m reactions 
A™- 1 -»• A™- 1 -»• ...A™- 1 (that are the cycle reac- 
tions without the limiting step) with corresponding 
constants from y™ -1 . 

If there exists an outgoing reaction A™ 1 — ► B in 
V m then we substitute it by the reaction A™^ 1 -> 
I? with the same constant, i.e. outgoing reactions 
A™ —>•... are reattached to the heads of the limiting 
steps (Fig. 6). Let us rearrange reactions from V m of 
the form B — > A™ . These reactions have prototypes 
in V m_1 (before the last gluing). We simply restore 
these reactions. If there exists a reaction A - n — > A" 1 
then we find the prototype in V™" 1 , A — > B, and 
substitute the reaction by A™ -1 — > £> with the same 
constant, as for Af -> A™ . 

After that step is performed, the vertices set is 
.A m_1 , but the reaction set differs from the reactions 
of the network V m_1 : the limiting steps of cycles are 
excluded and the outgoing reactions of glued cycles 
are included (reattached to the heads of the limiting 
steps). To make the next step, we select vertices of 
A™ 1 ^ 1 that are glued cycles from y TO_2 5 substitute 
these vertices by vertices of cycles, delete the lim- 
iting steps, attach outgoing reactions to the heads 
of the limiting steps, and for incoming reactions re- 
store their prototypes from y m_2 , and so on. 

After all, we restore all the glued cycles, and con- 
struct an acyclic reaction network on the set A. This 
acyclic network approximates relaxation of the net- 
work W. We call this system the dominant system 
of W and use notation dommod(W). 



In the simplest case, the dominant system is de- 
termined by the ordering of constants. But for suffi- 
ciently complex systems we need to introduce aux- 
iliary elementary reactions. They appear after cycle 
gluing and have monomial rate constants of the form 
k q = Yii kf, where q are integers, but not manda- 
tory positive. The dominant system depends on the 
place of these monomial values among the ordered 
constants. For systems with well separated constants 
we can also assume that each of these new constants 
will be well separated from other constants (Gorban 
& Radulescu (2008)). 

5.4. Example 

To demonstrate a possible branching of described 
algorithm for cycles surgery (gluing, restoring and 
cutting) with necessity of additional orderings, let 
us consider the following system: 

A 1 hA 2 hA s hA A hA 5 hA s , A 4 hA 2 , (31) 

(where the upper index marks the order of rate con- 
stants). The auxiliary discrete dynamical system for 
reaction network (31) is 

A 1 hA 2 hA 3 hA 4 h A 5 hA 3 . 

It has only one attractor, a cycle A 3 ^A^A^^A 3 . 
This cycle is not a sink for the whole network (31) be- 
cause reaction A 4 -^A 2 leads from that cycle. After 
gluing the cycle into a vertex A\ we get the new net- 
work Ai h-A 2 ^-A\ h-A 2 . The rate constant for the 
reaction A\^-A 2 is fc 23 = k 24 k 3 ^,/k^ 4 , where kij is 
the rate constant for the reaction Aj — > Ai in the 
initial network (/c 35 is the cycle limiting reaction). 
The new network coincides with its auxiliary system 
and has one cycle, A 2 ^A\^A 2 . This cycle is a sink, 
hence, we can start the back process of cycles restor- 
ing and cutting. One question arises immediately: 
which constant is smaller, k 32 or k 23 . The smallest 
of them is the limiting constant, and the answer de- 
pends on this choice. Let us consider two possibili- 
ties separately: (1) k 32 > k\ 3 and (2) fc 32 < k\ 3 . 

(1) Let as assume that k 32 > k 23 . The final auxil- 
iary system after gluing cycles is A\^A 2 ^)-A\^A 2 . 
Let us delete the limiting reaction A\^A 2 from the 
cycle. We get an acyclic system A\ : hA 2 ^rA\. The 
component A\ is the glued cycle A 3 ^hA 4 ^tAc,^-A 3 . 
Let us restore this cycle and delete the limiting 
reaction A 5 ^A 3 . We get the dominant system 
Ai^A^A^AihA^. Relaxation of this system 
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approximates relaxation of the initial network (31) 
under additional condition k 32 > k 23 - 

(2) Let as assume now that k 3 2 < ^23- The fi- 
nal auxiliary system after gluing cycles is the same, 
Ai-hA2^hAl^A2, but the limiting step in the cycle 
is different, A 2 ^A\. After cutting this step, we get 
acyclic system Ai^hA2^-A\, where the last reaction 
has rate constant k\ 3 ■ 

The component A\ is the glued cycle 

A 3 ^A A hA^A 3 . 

Let us restore this cycle and delete the limiting re- 
action A^^A 3 . The connection from glued cycle 
A\^fA2 with constant k 23 transforms into connec- 
tion A 5 ^hA 2 with the same constant k 23 - 
We get the dominant system: 

AihA 2 , A 3 ^tA A hA 5 hA 2 . 

The order of constants is now known: k 2 \ > £43 > 
fc 54 > k 23 , and we can substitute the sign "?" by 
"4": A 3 ^A 4 hA 5 ^A 2 . 

For both cases, k 32 > ^23 (^23 = ^24^35/^54) and 
^32 < ^23 ^ is easy to find the eigenvectors explicitly 
and to write the solution to the kinetic equations in 
explicit form. 

6. The Reversible Triangle of Reactions 

In this section, we illustrate the analysis of dom- 
inant systems on a simple example, the reversible 
triangle of reactions. 

A t o A 2 o A 3 A x (32) 

This triangle appeared in many works as an ideal 
object for a case study. Our favorite example is the 
work of Wei & Prater (1962). Now in our study the 
triangle (32) is not necessarily a closed system. We 
can assume that it is a subsystem of a larger sys- 
tem, and any reaction Ai — > Aj represents a reac- 
tion of the form . . . + Ai — »■ Aj + . . . , where unknown 
but slow components are substituted by dots. This 
means that there are no mandatory relations be- 
tween reaction rate constants, and six reaction rate 
constants are arbitrary nonnegative numbers. 

Let the reaction rate constant &21 for the reaction 
Ai — > A 2 be the largest. 

Let us describe all possible auxiliary dynamical 
systems for the triangle (32). For each vertex, we 
have to select the fastest outgoing reaction. For A\, 
it is always A\ — > A 2 , because of our choice of enu- 
meration (the higher scheme in Fig. 7). There exist 
two choices of the fastest outgoing reaction for two 




Fig. 7. Four possible auxiliary dynamical systems for the re- 
versible triangle of reactions with £121 > fcij for (i, j) ^ (2, 1): 
(a) fei2 > k Z 2, k 23 > k 13 ; (b) k 12 > k 32 , k 13 > k 23 ; (c) 
k:j 2 > k\ 2 , k 23 > £13; (d) £132 > fci2, £13 > k 23 . For each 
vertex the outgoing reaction with the largest rate constant 
is represented by the solid bold arrow, and other reactions 
are represented by the dashed arrows. The digraphs formed 
by solid bold arrows are the auxiliary discrete dynamical 
systems. Attractors of these systems are isolated in frames. 

other vertices and, therefore, only four versions of 
auxiliary dynamical systems for (32) (Fig. 7). Let 
us analyze in detail case (a). For the cases (b) and 
(c) the details of computations are similar. The ir- 
reversible cycle (d) is even simpler and was already 
discussed. 

6.1. Auxiliary System (a): Ai o A 2 <— A 3 ; 
ki2 > k 32 , k 23 > k 13 

6.1.1. Gluing Cycles 

The attractor is a cycle (with only two vertices) 
Ax ^ A 2 . This is not a sink, because two outgoing 
reactions exist: Ai — > A 3 and A 2 —¥ A 3 . They are 
relatively slow: k 3 \ <C /c 2 i and fc 32 <C k\2- The limit- 
ing step in this cycle is A 2 — > A\ with the rate con- 
stant ki2- We have to glue the cycle A\ <-> ^4 2 into 
one new component A\ and to add a new reaction 
A\ — > A 3 with the rate constant (see Fig. 5) 

&3J = max{/c 32 , k 31 k 12 /k 2 i} ■ (33) 

As a result, we get a new system, A\ «->• A 3 with 
reaction rate constants fc 31 (for A\ — > A 3 ) and initial 
fe3 (for A\ <— ^4 3 ). This cycle is a sink, because 
it has no outgoing reactions (the whole system is a 
trivial example of a sink) . 

6.1.2. Dominant System 

At the next step, we have to restore and cut the 
cycles. First cycle to cut is the result of cycle gluing, 
A\ f-> A 3 . It is necessary to delete the limiting step, 
i.e. the reaction with the smallest rate constant. If 
^31 > &23, then we get A\ ->A 3 . If, inverse, k 23 > 
k 31 , then we obtain A\ <— A 3 . 
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Fig. 8. Dominant systems for case (a) (defined in Fig. 7) 

After that, we have to restore and cut the cycle 
which was glued into the vertex A \. This is the two- 
vertices cycle A\ «-> A 2 . The limiting step for this 
cycle is A\ <— A 2 , because k 2 \ ~> k\ 2 - If £31 > k 23 , 
then following the rule visualized by Fig. 6, we get 
the dominant system A\ — > A 2 — > A3 with reaction 
rate constants k 2 \ for A\ — > A 2 and k 31 for A 2 — > 
A3. If k 2 3 > k 31 then we obtain Ai ^ A 2 <— A 3 with 
reaction rate constants fei for A\ — > A 2 and fc 23 for 
A 2 ^— A3. All the procedure is illustrated by Fig. 8. 



6.1.3. Eigenvalues and Eigenvectors 

The eigenvalues and the corresponding eigenvec- 
tors for dominant systems in case (a) are represented 
below in zero-one asymptotic. 

(i) k^ > k 23 , 

the dominant system A\ — > A 2 — > As, 



A o = 0, r°» (0,0,1), i° = (1,1,1) 

Ai w -k 21 , r 1 « (1,-1,0), Z 1 w (1,0,0) 

A2«-*Ji, r 2 « (0,1,-1), Z 2 « (1,1,0) 

(34) 

(ii) fc 23 > fcji, 

the dominant system Ai — > A 2 <— A3, 

A =0, r°» (0,1,0), i° = (1,1,1); 



Ai 



-fc 21 , r 1 ** (1,-1,0), Z 1 - (1,0,0); 



A 2 «-fc 23 , r 2 « (0,-1,1), Z 2 « (0,0,1). 

(35) 

Here, the value of fc^ is given by formula (33). 

Analysis of examples provided us by an impor- 
tant conclusion: the number of different dominant 
systems in examples was less than the number of 
all possible orderings. For many pairs of constants 
kij, ki r it is not important which of them is larger. 
There is no need to consider all orderings of mono- 
mials. We have to consider only those inequalities 
between constants and monomials that appear in 
the construction of the dominant systems. 



7. Corrections to Dominant Dynamics 

The hierarchy of systems W, W 1 , W 2 , ... can be 
used for multigrid correction of the dominant dy- 
namics. The simple example of multigrid approach 
gives the algorithm of steady state approximation 
(Gorban & Radulescu (2008)). For this purpose, on 
the way up (cycle restoration and cutting, Sec. 5.3.2) 
we calculate distribution in restoring cycles with 
higher accuracy, by exact formula (13), or in linear 
approximation (15) instead of the simplest zero-one 
asymptotic (16). Essentially, the way up remains the 
same. 

After termination of the gluing process, we can 
find all steady state distributions by restoring cy- 
cles in the auxiliary reaction network V m . Let 
A™ , A™ , ... be fixed points of $ m . The set of steady 
states for V m is the set of all distributions on the 
set of fixed points {A™, A™ , ...}. 

Let us take one of the basis distributions, = 1, 
other a = on V m . If the vertex A^ is a glued cycle, 
then we substitute them by all the vertices of this cy- 
cle. Redistribute the concentration cjl between the 
vertices of the corresponding cycle by the rule (13) 
(or by an approximation). As a result, we get a set 
of vertices and a distribution on this set of vertices. 
If among these vertices there are glued cycles, then 
we repeat the procedure of cycle restoration. Termi- 
nate when there is no glued cycles in the support of 
the distribution. 

The resulting distribution is the approximation to 
a steady state of W, and the basis of steady states 
for W can be approximated by this method. 

For example, for the system Fig. 8 we have, first 
of all, to compute the stationary distribution in the 
cycle A\ o A3, c\ and c 3 . On the base of the general 
formula for a simple cycle (13) we obtain: 



w = 



4- _L_ 



k 1 
K 31 



w 

k 23 



Tj~ ' C3 = T, • ( 36 ) 



After that, we have to restore the cycle glued into 
A\. This means to calculate the concentrations of 



A\ and A 2 with normalization c\ - 
(13) gives: 



c 2 = c\. Formula 



w 1 = j I 1 , Cl = — , c 2 = — . (37) 

For eigenvectors, there appear two operations of 
corrections: (i) correction for an acyclic network 
without branching (43), (45), and (ii) corrections for 
a cycle with relatively slow outgoing reactions (49). 
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These corrections are by-products of the accuracy 
estimates given in Appendix. 

8. Conclusion 

Now, the idea of limiting step is developed to the 
asymptotology of multiscale reaction networks. We 
found the main terms of eigenvectors and eigen- 
values asymptotic on logarithmic straight lines 
lnfcjj = when £ — > oo. These main terms could 
be represented by acyclic dominant system which 
is a piecewise constant function of the direction 
vectors (Oij). This theory gives the analogue of 
the Vishik & Ljusternik (1960) theory for chemical 
reaction networks. We demonstrated also how to 
construct the accuracy estimates and the first order 
corrections to eigenvalues and eigenvectors. 

There are several ways of using the developed the- 
ory and algorithms: 

- For direct computation of steady states and relax- 
ation dynamics; this may be useful for complex 
systems because of the simplicity of the algorithm 
and resulting formulas and because often we do 
not know the rate constants for complex networks, 
and kinetics that is ruled by orderings rather than 
by exact values of rate constants may be very use- 
ful in practically frequent situation when the val- 
ues of the various reaction constants are unknown 
or poorly known; 

- For planning experiments and mining the exper- 
imental data - the observable kinetics is more 
sensitive to reactions from the dominant net- 
work, and much less sensitive to other reactions, 
the relaxation spectrum of the dominant network 
is explicitly connected with the correspondent 
reaction rate constants, and the eigenvectors 
( "modes" ) are sensitive to the constant ordering, 
but not to exact values; 

- The steady states and dynamics of the dominant 
system could serve as a robust first approximation 
in perturbation theory or as a preconditioning in 
numerical methods. 

The next step should be development of asymp- 
totic estimates for networks with modular struc- 
ture and time separations between modules, not be- 
tween individual reactions. But now it seems that 
the most important further development should be 
the asymptotology of nonlinear reaction networks. 
For multiscale nonlinear reaction networks the ex- 
pected dynamical behaviour is to be approximated 
by the system of dominant networks. These net- 



works may change in time (this is the significant dif- 
ference from the linear case) but remain relatively 
simple. 
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Appendix: Mathematical Backgrounds of Ac- 
curacy Estimation 

Estimates for Perturbed Acyclic Networks 

The famous Gerschgorin theorem (Marcus & 
Mine (1992), Varga (2004)) gives estimates of 
eigenvalues. We need also estimates of eigenvec- 
tors. Below A = (ciij) is a complex n x n matrix, 
Qi = J2j jjii \ a ji\ (sums of non-diagonal elements 
in columns). 

Gerschgorin theorem (Marcus & Mine (1992), 
p. 146): The characteristic roots of A lie in the closed 
region G Q of the z-planc 

G Q = \jG? (G? = {z\\z-a ii \<Q i }. (38) 

i 

Areas Gf are the Gerschgorin discs. (The same esti- 
mate are valid for sums in rows, Pj. Here and below 
we don't duplicate the estimates.) 

Gerschgorin disks Gf (i = 1, . . . n) are isolated, 
if Gf n Gf = for i ^ j. If disks Gf (i = 1, . . . n) 
are isolated, then the spectrum of A is simple, and 
each Gerschgorin disk Gf contains one and only one 
eigenvalue of A (Marcus & Mine (1992), p. 147). 

We assume that Gerschgorin disks Gf (i = 
l,...n) are isolated: for all i, j (i ^ j) 

\a ii -a jj \> Qi + Qj. (39) 

Let us introduce the following notations: 

Qi _ \ai ' 

i i — 

a,;,; 



max Si 



Xij > X : 



max Xij, 



mm 



^33 I 



g = mm.g,. 



(40) 

Usually, we consider Ei and Xij as sufficiently small 
numbers. In contrary, the diagonal gap g should not 
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be small, (this is the gap condition). For example, if 
for any two diagonal elements an, ajj either an 3> 
ajj or an <C ajj, then > 1 for all i. 

Let Xi e Gf be the eigenvalue of A (|A, — on| < 
Qi). Let us estimate the corresponding right eigen- 
vector r-W . We take r\ = 1 and for j 7^ i introduce 



a (n — l)-dimensional vector x 1 : x^ 
(i ^ j). For x l we get equation 

(1 = -5* 



r!.(a 



(41) 



where a 1 is a vector of the non-diagonal elements 
of the ith column of A (a* = a^ , j 7^ i), and the 
(n— l)x(n — 1) matrix B % has matrix elements 
(3,1 ^i) 

V ; " °> l -(i^j) (42) 



(0 



< 



Due to the Gerschgorin estimate, \b 
From Eq. (41) we obtain: 

x i = -a i -BW(l-B®)- 1 a i . (43) 

From this definition and simple estimates in I 1 norm, 
we get the following estimate of eigenvectors. 

Theorem 2. Let the Gerschgoring disks be iso- 
lated, and the diagonal gap be big enough: g > ne. 
Then for the zth eigenvector of A the following uni- 
form estimate holds: 

-2 



x. 



ne 



(j ? 1, r\ = 1). □ (44) 

g(g - ne) 

So, if the matrix A is diagonally dominant and the 
diagonal gap g is big enough, then the eigenvectors 
are proven to be close to the standard basis vectors 
with explicit evaluation of accuracy. 

The first correction to eigenvectors is also given by 
Eq. (43). If for the iteration we use the Gerschgorin 
estimates for eigenvalue Aj ~ an , then we can write 
in the next approximation for eigenvectors (r- = 



*3i 



(45) 



where B^j is the non-diagonal part of : it has 
the same non-diagonal elements and zeros on diago- 
nal. There exists plenty of further simplifications for 
this iteration formula. For example, one can leave 
just the first term, that gives the first order approx- 
imation in the power of £ (x < e). 

To apply these estimates to an acyclic network 
supplemented by additional reactions, we have to 
use the eigenbasis of this acyclic network (Sec. 4). 



Direct use of this theorem and estimates for a kinetic 
matrix K in the standard basis is impossible, the 
diagonal dominance in this coordinate system is not 
large, and sums of elements in columns are zero. To 
apply this theorem we need two lemmas. 

Let W be a reaction network without branching 
(a finite dynamical system) with n vertices. Then 
the number of reactions in W is n — /, where / is 
the number of fixed points (the vertices without out- 
going reactions) . Let Y be the set of stoichiometric 
vectors for W. 

Lemma 1. Y forms a basis in the subspace 
{c I J2i °i — 0} ^ an d on ly if the reaction network 
W is acyclic and connected (has only one fixed 
point). □ 

Let us consider a general reaction network on the 
set Ai,...A n . For stoichiometric vector of reaction 
Ai —> Ai we use notation 7^. Assume that the auxil- 
iary dynamical system i 1— > (f>(i) for a given reaction 
network is acyclic and has only one attractor, a fixed 
point. For this auxiliary network, we use notation: 
«i = kji for the only reaction Ai — > Aj, or Ki = 0. 

For every reaction of the initial network, Ai —¥ A{, 
a linear operators Qu can be defined by its action 
on the basis vectors, "y<j,(i)i- 

Qil{l<t,{i) i) = Hi, Qu(l4,( P ) P ) = for p ^ i. (46) 
Lemma 2. The kinetic equation for the whole 

reaction network (9) could be transformed to the 

form 



1+ —Q]i\^Zl<t>(i)i^ (47) 



3,1 d^M.3)) 



= 1+ E 



Kc, 



3,1 



where K is kinetic matrix of the kinetic equation for 
the auxiliary network. □ 

By construction of auxiliary dynamical system, 
kn < Ki if I 7^ 4>(i), and for reaction networks with 
well separated constants ku <C «». Notice also that 
the matrix Qji docs not depend on rate constants 
values. 

For matrix K we have the eigenbasis in explicit 
form. Let us represent system (47) in this eigenbasis 
of K. Any matrix B in this eigenbasis has the form 
B = (bij), bij = l l Br 1 = Y, qs l q bqsr J s , where (b qs ) is 
matrix B in the initial basis, V and r J are left and 
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right eigenvectors of K (27), (28). In eigenbasis of K 
the estimates of eigenvalues and estimates of eigen- 
vectors are much more efficient than in original coor- 
dinates: the system is strongly diagonally dominant. 
Transformation to this basis is an effective precondi- 
tioning for the perturbation theory that uses auxil- 
iary kinetics as a first approximation to the kinetics 
of the whole system. 

Estimates for Perturbed Ergodic Systems 

Let us consider a strongly connected network 
with kinetic matrix K. The corresponding kinet- 
ics is ergodic and there exists unique normalized 
steady state c* > 0, c* = 1. For each i we define 
Ki = kji. The number — Ki is the nth diagonal 
element of unperturbed kinetic matrix K. 

Let this network be perturbed by outgoing reac- 
tions Ai — > 0. The perturbation has the "loss form" : 
the perturbed matrix is K — diag(£i m ) , perturbation 
of each diagonal element is relatively small (diag is 
the diagonal matrix) . 

The perturbations e relatively small with 

respect to k i , but not obligatory small with respect 
to other rate constants. 

First, we do not assume anything about value of 
£i > and make the following transformation. For 
an arbitrary normalized vector r (r^ > 0, r, = I) 
we add to the network reactions Ai — > Aj with reac- 
tion rates qji = rjSiHi. We use Q(r) for the kinetic 
matrix of this additional network. Simple algebra 
gives 

Q(r) + diag(e 4 K l ) = [e^r, s 2 K2r, ...e„K„r] 
= r(£iKi,e 2 «2, —e n n n ). 

Here, in the right hand side we have a matrix, all 
columns of which are proportional to the vector r, 
this is a product of r on the vector-raw of coefficients. 
We represent the perturbed matrix in the form K — 
diag(eiKi) = K + Q(r) - (Q(r) + diag^Kj)). 

Theorem 3. There exists such normalized posi- 
tive r* that (K + Q(r*))r* = 0. This r* is an eigen- 
vector of the perturbed network with the eigenvalue 
A = J2i r i £ i K i, an d, at the same time, it is a steady- 
state for the network with kinetic matrix K + Q(r*). 

To prove existence it is sufficient to mention, that 
for any r the network with kinetic matrix K + Q(r) 
has unique positive normalized steady state c*(r), 
which depends continuously on r. The map r M> 
c* (r) has a fixed point r* (the Brouwer fixed point 
theorem). □ 



This representation allows us to produce useful 
estimates, for example, when the unperturbed sys- 
tem is a cycle, we find \r* — c* < 3e|c*| under con- 
dition e < 0.25, where e = J2 £ i- Formula for the 
first correction gives (r* = c* + <5r,, w = he*): 

i 

( 49 ) 

v = — > 2(ec 4 - Si), 
n ' 

i=l 

For more complex networks, the explicit formulas 
for corrections could be produced on the base of the 
network graphs, similar to the steady-state formu- 
las, presented, for example, by Yablonskii, Bykov, 
Gorban, & Elokhin (1991). 

So, the asymptotic analysis gives good approxi- 
mation of eigenvectors and eigenvalues for kinetic 
matrix. The condition number is big (unbounded) 
but these estimates work even better when the con- 
stants become more separated. Nevertheless, some 
caution is needed: the error is proven to be small, 
but the residuals (the values \\Kr — Ar|| for approx- 
imations of r and A) may be not small (Gorban & 
Radulescu (2008)). 
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